

<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
  <meta charset="utf-8">
  
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  
  <title>Tutorial of how to train a GAP model from scratch &mdash; quippy 7ea04c5+ documentation</title>
  

  
  
    <link rel="shortcut icon" href="../_static/favicon.ico"/>
  
  
  

  
  <script type="text/javascript" src="../_static/js/modernizr.min.js"></script>
  
    
      <script type="text/javascript" id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
        <script type="text/javascript" src="../_static/jquery.js"></script>
        <script type="text/javascript" src="../_static/underscore.js"></script>
        <script type="text/javascript" src="../_static/doctools.js"></script>
        <script type="text/javascript" src="../_static/language_data.js"></script>
        <script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
        <script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"processClass": "math|output_area", "processEscapes": true, "ignoreClass": "document", "inlineMath": [["$", "$"], ["\\(", "\\)"]]}})</script>
    
    <script type="text/javascript" src="../_static/js/theme.js"></script>

    

  
  <link rel="stylesheet" href="../_static/css/theme.css" type="text/css" />
  <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
  <link rel="stylesheet" href="../_static/theme_overrides.css" type="text/css" />
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" /> 
</head>

<body class="wy-body-for-nav">

   
  <div class="wy-grid-for-nav">
    
    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
      <div class="wy-side-scroll">
        <div class="wy-side-nav-search" >
          

          
            <a href="../index.html" class="icon icon-home"> quippy
          

          
            
            <img src="../_static/hybrid.png" class="logo" alt="Logo"/>
          
          </a>

          
            
            
              <div class="version">
                7ea04c5+
              </div>
            
          

          
<div role="search">
  <form id="rtd-search-form" class="wy-form" action="../search.html" method="get">
    <input type="text" name="q" placeholder="Search docs" />
    <input type="hidden" name="check_keywords" value="yes" />
    <input type="hidden" name="area" value="default" />
  </form>
</div>

          
        </div>

        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
          
    
            
            
              
            
            
              <ul>
<li class="toctree-l1"><a class="reference internal" href="../intro.html">Introduction to <code class="docutils literal notranslate"><span class="pre">quippy</span></code></a></li>
<li class="toctree-l1"><a class="reference internal" href="../install.html">Installation of QUIP and quippy</a></li>
<li class="toctree-l1"><a class="reference internal" href="index.html">Tutorials</a></li>
<li class="toctree-l1"><a class="reference internal" href="../quippy.html">Quippy library reference</a></li>
<li class="toctree-l1"><a class="reference internal" href="../visualisation.html">Visualisation Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../fortran_wrapper.html">Appendix: wrapping Fortran 90 Code</a></li>
<li class="toctree-l1"><a class="reference internal" href="../teach_sparse.html">Command line options of the teach_sparse main program</a></li>
</ul>

            
          
    <p></p>
    <ul>
        <li><a href="../genindex.html">Index</a></li>
        <li><a href="../py-modindex.html">Module Index</a></li>
        <li><a href="../search.html">Search Page</a></li>
        <li><a href="../_modules/index.html">Module Source Code</a></li>
    </ul>
  
        </div>
      </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">

      
      <nav class="wy-nav-top" aria-label="top navigation">
        
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="../index.html">quippy</a>
        
      </nav>


      <div class="wy-nav-content">
        
        <div class="rst-content">
        
          















<div role="navigation" aria-label="breadcrumbs navigation">

  <ul class="wy-breadcrumbs">
    
      <li><a href="../index.html">Docs</a> &raquo;</li>
        
      <li>Tutorial of how to train a GAP model from scratch</li>
    
    
      <li class="wy-breadcrumbs-aside">
        
            
            <a href="../_sources/Tutorials/gap_fitting_tutorial.ipynb.txt" rel="nofollow"> View page source</a>
          
        
      </li>
    
  </ul>

  
  <hr/>
</div>
          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
           <div itemprop="articleBody">
            
  
<style>
/* CSS for nbsphinx extension */

/* remove conflicting styling from Sphinx themes */
div.nbinput,
div.nbinput div.prompt,
div.nbinput div.input_area,
div.nbinput div[class*=highlight],
div.nbinput div[class*=highlight] pre,
div.nboutput,
div.nbinput div.prompt,
div.nbinput div.output_area,
div.nboutput div[class*=highlight],
div.nboutput div[class*=highlight] pre {
    background: none;
    border: none;
    padding: 0 0;
    margin: 0;
    box-shadow: none;
}

/* avoid gaps between output lines */
div.nboutput div[class*=highlight] pre {
    line-height: normal;
}

/* input/output containers */
div.nbinput,
div.nboutput {
    display: -webkit-flex;
    display: flex;
    align-items: flex-start;
    margin: 0;
    width: 100%;
}
@media (max-width: 540px) {
    div.nbinput,
    div.nboutput {
        flex-direction: column;
    }
}

/* input container */
div.nbinput {
    padding-top: 5px;
}

/* last container */
div.nblast {
    padding-bottom: 5px;
}

/* input prompt */
div.nbinput div.prompt pre {
    color: #307FC1;
}

/* output prompt */
div.nboutput div.prompt pre {
    color: #BF5B3D;
}

/* all prompts */
div.nbinput div.prompt,
div.nboutput div.prompt {
    min-width: 5ex;
    padding-top: 0.4em;
    padding-right: 0.4em;
    text-align: right;
    flex: 0;
}
@media (max-width: 540px) {
    div.nbinput div.prompt,
    div.nboutput div.prompt {
        text-align: left;
        padding: 0.4em;
    }
    div.nboutput div.prompt.empty {
        padding: 0;
    }
}

/* disable scrollbars on prompts */
div.nbinput div.prompt pre,
div.nboutput div.prompt pre {
    overflow: hidden;
}

/* input/output area */
div.nbinput div.input_area,
div.nboutput div.output_area {
    padding: 0.4em;
    -webkit-flex: 1;
    flex: 1;
    overflow: auto;
}
@media (max-width: 540px) {
    div.nbinput div.input_area,
    div.nboutput div.output_area {
        width: 100%;
    }
}

/* input area */
div.nbinput div.input_area {
    border: 1px solid #e0e0e0;
    border-radius: 2px;
    background: #f5f5f5;
}

/* override MathJax center alignment in output cells */
div.nboutput div[class*=MathJax] {
    text-align: left !important;
}

/* override sphinx.ext.imgmath center alignment in output cells */
div.nboutput div.math p {
    text-align: left;
}

/* standard error */
div.nboutput div.output_area.stderr {
    background: #fdd;
}

/* ANSI colors */
.ansi-black-fg { color: #3E424D; }
.ansi-black-bg { background-color: #3E424D; }
.ansi-black-intense-fg { color: #282C36; }
.ansi-black-intense-bg { background-color: #282C36; }
.ansi-red-fg { color: #E75C58; }
.ansi-red-bg { background-color: #E75C58; }
.ansi-red-intense-fg { color: #B22B31; }
.ansi-red-intense-bg { background-color: #B22B31; }
.ansi-green-fg { color: #00A250; }
.ansi-green-bg { background-color: #00A250; }
.ansi-green-intense-fg { color: #007427; }
.ansi-green-intense-bg { background-color: #007427; }
.ansi-yellow-fg { color: #DDB62B; }
.ansi-yellow-bg { background-color: #DDB62B; }
.ansi-yellow-intense-fg { color: #B27D12; }
.ansi-yellow-intense-bg { background-color: #B27D12; }
.ansi-blue-fg { color: #208FFB; }
.ansi-blue-bg { background-color: #208FFB; }
.ansi-blue-intense-fg { color: #0065CA; }
.ansi-blue-intense-bg { background-color: #0065CA; }
.ansi-magenta-fg { color: #D160C4; }
.ansi-magenta-bg { background-color: #D160C4; }
.ansi-magenta-intense-fg { color: #A03196; }
.ansi-magenta-intense-bg { background-color: #A03196; }
.ansi-cyan-fg { color: #60C6C8; }
.ansi-cyan-bg { background-color: #60C6C8; }
.ansi-cyan-intense-fg { color: #258F8F; }
.ansi-cyan-intense-bg { background-color: #258F8F; }
.ansi-white-fg { color: #C5C1B4; }
.ansi-white-bg { background-color: #C5C1B4; }
.ansi-white-intense-fg { color: #A1A6B2; }
.ansi-white-intense-bg { background-color: #A1A6B2; }

.ansi-default-inverse-fg { color: #FFFFFF; }
.ansi-default-inverse-bg { background-color: #000000; }

.ansi-bold { font-weight: bold; }
.ansi-underline { text-decoration: underline; }

/* Some additional styling taken form the Jupyter notebook CSS */
div.rendered_html table {
  border: none;
  border-collapse: collapse;
  border-spacing: 0;
  color: black;
  font-size: 12px;
  table-layout: fixed;
}
div.rendered_html thead {
  border-bottom: 1px solid black;
  vertical-align: bottom;
}
div.rendered_html tr,
div.rendered_html th,
div.rendered_html td {
  text-align: right;
  vertical-align: middle;
  padding: 0.5em 0.5em;
  line-height: normal;
  white-space: normal;
  max-width: none;
  border: none;
}
div.rendered_html th {
  font-weight: bold;
}
div.rendered_html tbody tr:nth-child(odd) {
  background: #f5f5f5;
}
div.rendered_html tbody tr:hover {
  background: rgba(66, 165, 245, 0.2);
}

/* CSS overrides for sphinx_rtd_theme */

/* 24px margin */
.nbinput.nblast,
.nboutput.nblast {
    margin-bottom: 19px;  /* padding has already 5px */
}

/* ... except between code cells! */
.nblast + .nbinput {
    margin-top: -19px;
}

.admonition > p:before {
    margin-right: 4px;  /* make room for the exclamation icon */
}

/* Fix math alignment, see https://github.com/rtfd/sphinx_rtd_theme/pull/686 */
.math {
    text-align: unset;
}
</style>
<div class="section" id="Tutorial-of-how-to-train-a-GAP-model-from-scratch">
<h1>Tutorial of how to train a GAP model from scratch<a class="headerlink" href="#Tutorial-of-how-to-train-a-GAP-model-from-scratch" title="Permalink to this headline">¶</a></h1>
<div class="section" id="steps">
<h2>steps<a class="headerlink" href="#steps" title="Permalink to this headline">¶</a></h2>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="mf">1.</span> <span class="n">generate</span> <span class="n">a</span> <span class="n">small</span> <span class="n">dataset</span> <span class="n">of</span> <span class="n">water</span> <span class="n">structures</span>
    <span class="o">-</span> <span class="n">use</span> <span class="n">CP2K</span> <span class="k">if</span> <span class="n">you</span> <span class="n">havea</span> <span class="n">access</span> <span class="n">to</span> <span class="n">it</span>
    <span class="o">-</span> <span class="n">otherwise</span><span class="p">:</span> <span class="n">use</span> <span class="nb">any</span> <span class="n">simple</span> <span class="n">potential</span> <span class="n">implemented</span> <span class="ow">in</span> <span class="n">ASE</span><span class="p">,</span> <span class="n">just</span> <span class="k">for</span> <span class="n">trying</span> <span class="n">this</span> <span class="n">out</span> <span class="n">I</span> <span class="n">have</span> <span class="n">used</span> <span class="n">EMT</span> <span class="n">here</span>
<span class="mf">1.</span> <span class="n">generate</span> <span class="n">e0</span> <span class="n">values</span>
<span class="mf">1.</span> <span class="n">separate</span> <span class="n">a</span> <span class="n">training</span> <span class="ow">and</span> <span class="n">a</span> <span class="n">validation</span> <span class="n">dataset</span>
<span class="mf">1.</span> <span class="o">**</span><span class="n">train</span> <span class="n">the</span> <span class="n">model</span><span class="o">**</span>
<span class="mf">1.</span> <span class="n">look</span> <span class="n">at</span> <span class="n">the</span> <span class="n">outcome</span> <span class="n">of</span> <span class="n">the</span> <span class="n">model</span>
</pre></div>
</div>
</div>
<div class="section" id="here-we-will-fit-twice,-to-see-the-difference-between-a-2b-only-and-a-2b+3b-fit">
<h2>here we will fit twice, to see the difference between a 2b-only and a 2b+3b fit<a class="headerlink" href="#here-we-will-fit-twice,-to-see-the-difference-between-a-2b-only-and-a-2b+3b-fit" title="Permalink to this headline">¶</a></h2>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[15]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="c1"># general imports</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">from</span> <span class="nn">copy</span> <span class="k">import</span> <span class="n">deepcopy</span> <span class="k">as</span> <span class="n">cp</span>

<span class="c1"># ase imports</span>
<span class="kn">import</span> <span class="nn">ase.io</span>
<span class="kn">from</span> <span class="nn">ase</span> <span class="k">import</span> <span class="n">Atoms</span><span class="p">,</span> <span class="n">Atom</span>
<span class="kn">from</span> <span class="nn">ase</span> <span class="k">import</span> <span class="n">units</span>
<span class="kn">from</span> <span class="nn">ase.build</span> <span class="k">import</span> <span class="n">molecule</span>
<span class="c1"># for MD</span>
<span class="kn">from</span> <span class="nn">ase.md.langevin</span> <span class="k">import</span> <span class="n">Langevin</span>
<span class="kn">from</span> <span class="nn">ase.io.trajectory</span> <span class="k">import</span> <span class="n">Trajectory</span>
</pre></div>
</div>
</div>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[2]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="c1"># helper functions</span>
<span class="k">def</span> <span class="nf">make_water</span><span class="p">(</span><span class="n">density</span><span class="p">,</span> <span class="n">super_cell</span><span class="o">=</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">]):</span>
    <span class="sd">&quot;&quot;&quot; Geenrates a supercell of water molecules with a desired density.</span>
<span class="sd">        Density in g/cm^3!!!&quot;&quot;&quot;</span>
    <span class="n">h2o</span> <span class="o">=</span> <span class="n">molecule</span><span class="p">(</span><span class="s1">&#39;H2O&#39;</span><span class="p">)</span>
    <span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cbrt</span><span class="p">((</span><span class="nb">sum</span><span class="p">(</span><span class="n">h2o</span><span class="o">.</span><span class="n">get_masses</span><span class="p">())</span> <span class="o">*</span> <span class="n">units</span><span class="o">.</span><span class="n">m</span> <span class="o">**</span> <span class="mi">3</span> <span class="o">*</span> <span class="mf">1E-6</span> <span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">density</span> <span class="o">*</span> <span class="n">units</span><span class="o">.</span><span class="n">mol</span><span class="p">))</span>
    <span class="n">h2o</span><span class="o">.</span><span class="n">set_cell</span><span class="p">((</span><span class="n">a</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">a</span><span class="p">))</span>
    <span class="n">h2o</span><span class="o">.</span><span class="n">set_pbc</span><span class="p">((</span><span class="kc">True</span><span class="p">,</span> <span class="kc">True</span><span class="p">,</span> <span class="kc">True</span><span class="p">))</span>
    <span class="c1">#return cp(h2o.repeat(super_cell))</span>
    <span class="k">return</span> <span class="n">h2o</span><span class="o">.</span><span class="n">repeat</span><span class="p">(</span><span class="n">super_cell</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">rms_dict</span><span class="p">(</span><span class="n">x_ref</span><span class="p">,</span> <span class="n">x_pred</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot; Takes two datasets of the same shape and returns a dictionary containing RMS error data&quot;&quot;&quot;</span>

    <span class="n">x_ref</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">x_ref</span><span class="p">)</span>
    <span class="n">x_pred</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">x_pred</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">shape</span><span class="p">(</span><span class="n">x_pred</span><span class="p">)</span> <span class="o">!=</span> <span class="n">np</span><span class="o">.</span><span class="n">shape</span><span class="p">(</span><span class="n">x_ref</span><span class="p">):</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;WARNING: not matching shapes in rms&#39;</span><span class="p">)</span>

    <span class="n">error_2</span> <span class="o">=</span> <span class="p">(</span><span class="n">x_ref</span> <span class="o">-</span> <span class="n">x_pred</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span>

    <span class="n">average</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="n">error_2</span><span class="p">))</span>
    <span class="n">std_</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">var</span><span class="p">(</span><span class="n">error_2</span><span class="p">))</span>

    <span class="k">return</span> <span class="p">{</span><span class="s1">&#39;rmse&#39;</span><span class="p">:</span> <span class="n">average</span><span class="p">,</span> <span class="s1">&#39;std&#39;</span><span class="p">:</span> <span class="n">std_</span><span class="p">}</span>
</pre></div>
</div>
</div>
</div>
<div class="section" id="generating-data-only-with-ASE,-using-the-EMT-calculator">
<h2>generating data only with ASE, using the EMT calculator<a class="headerlink" href="#generating-data-only-with-ASE,-using-the-EMT-calculator" title="Permalink to this headline">¶</a></h2>
<p>This is only for the demonstration of how to do it, this run is will be done very fast. There is no practical use of the data beyond lerning the use teach_sparse, quip, etc. with it.</p>
<div class="nbinput docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[3]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="c1"># Running MD with ASE&#39;s EMT</span>

<span class="kn">from</span> <span class="nn">ase.calculators.emt</span> <span class="k">import</span> <span class="n">EMT</span>
<span class="n">calc</span> <span class="o">=</span> <span class="n">EMT</span><span class="p">()</span>

<span class="n">T</span> <span class="o">=</span> <span class="mi">150</span>  <span class="c1"># Kelvin</span>

<span class="c1"># Set up a grid of water</span>
<span class="n">water</span> <span class="o">=</span> <span class="n">make_water</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">])</span>
<span class="n">water</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">calc</span><span class="p">)</span>

<span class="c1"># We want to run MD using the Langevin algorithm</span>
<span class="c1"># with a time step of 1 fs, the temperature T and the friction</span>
<span class="c1"># coefficient to 0.002 atomic units.</span>
<span class="n">dyn</span> <span class="o">=</span> <span class="n">Langevin</span><span class="p">(</span><span class="n">water</span><span class="p">,</span> <span class="mi">1</span> <span class="o">*</span> <span class="n">units</span><span class="o">.</span><span class="n">fs</span><span class="p">,</span> <span class="n">T</span> <span class="o">*</span> <span class="n">units</span><span class="o">.</span><span class="n">kB</span><span class="p">,</span> <span class="mf">0.0002</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">printenergy</span><span class="p">(</span><span class="n">a</span><span class="o">=</span><span class="n">water</span><span class="p">):</span>  <span class="c1"># store a reference to atoms in the definition.</span>
    <span class="sd">&quot;&quot;&quot;Function to print the potential, kinetic and total energy.&quot;&quot;&quot;</span>
    <span class="n">epot</span> <span class="o">=</span> <span class="n">a</span><span class="o">.</span><span class="n">get_potential_energy</span><span class="p">()</span> <span class="o">/</span> <span class="nb">len</span><span class="p">(</span><span class="n">a</span><span class="p">)</span>
    <span class="n">ekin</span> <span class="o">=</span> <span class="n">a</span><span class="o">.</span><span class="n">get_kinetic_energy</span><span class="p">()</span> <span class="o">/</span> <span class="nb">len</span><span class="p">(</span><span class="n">a</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Energy per atom: Epot = </span><span class="si">%.3f</span><span class="s1">eV  Ekin = </span><span class="si">%.3f</span><span class="s1">eV (T=</span><span class="si">%3.0f</span><span class="s1">K)  &#39;</span>
          <span class="s1">&#39;Etot = </span><span class="si">%.3f</span><span class="s1">eV&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">epot</span><span class="p">,</span> <span class="n">ekin</span><span class="p">,</span> <span class="n">ekin</span> <span class="o">/</span> <span class="p">(</span><span class="mf">1.5</span> <span class="o">*</span> <span class="n">units</span><span class="o">.</span><span class="n">kB</span><span class="p">),</span> <span class="n">epot</span> <span class="o">+</span> <span class="n">ekin</span><span class="p">))</span>

<span class="n">dyn</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">printenergy</span><span class="p">,</span> <span class="n">interval</span><span class="o">=</span><span class="mi">5</span><span class="p">)</span>

<span class="c1"># We also want to save the positions of all atoms after every 5th time step.</span>
<span class="n">traj</span> <span class="o">=</span> <span class="n">Trajectory</span><span class="p">(</span><span class="s1">&#39;dyn_emt.traj&#39;</span><span class="p">,</span> <span class="s1">&#39;w&#39;</span><span class="p">,</span> <span class="n">water</span><span class="p">)</span>
<span class="n">dyn</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">traj</span><span class="o">.</span><span class="n">write</span><span class="p">,</span> <span class="n">interval</span><span class="o">=</span><span class="mi">5</span><span class="p">)</span>

<span class="c1"># Now run the dynamics</span>
<span class="n">printenergy</span><span class="p">(</span><span class="n">water</span><span class="p">)</span>
<span class="n">dyn</span><span class="o">.</span><span class="n">run</span><span class="p">(</span><span class="mi">600</span><span class="p">)</span>   <span class="c1"># CHANGE THIS IF YOU WANT LONGER/SHORTER RUN</span>
</pre></div>
</div>
</div>
<div class="nboutput nblast docutils container">
<div class="prompt empty docutils container">
</div>
<div class="output_area docutils container">
<div class="highlight"><pre>
Energy per atom: Epot = 0.885eV  Ekin = 0.000eV (T=  0K)  Etot = 0.885eV
Energy per atom: Epot = 0.820eV  Ekin = 0.053eV (T=413K)  Etot = 0.874eV
Energy per atom: Epot = 0.660eV  Ekin = 0.208eV (T=1611K)  Etot = 0.868eV
Energy per atom: Epot = 0.632eV  Ekin = 0.224eV (T=1736K)  Etot = 0.857eV
Energy per atom: Epot = 0.869eV  Ekin = 0.005eV (T= 39K)  Etot = 0.874eV
Energy per atom: Epot = 0.781eV  Ekin = 0.088eV (T=682K)  Etot = 0.869eV
Energy per atom: Epot = 0.756eV  Ekin = 0.117eV (T=903K)  Etot = 0.872eV
Energy per atom: Epot = 0.706eV  Ekin = 0.164eV (T=1270K)  Etot = 0.870eV
Energy per atom: Epot = 0.703eV  Ekin = 0.158eV (T=1223K)  Etot = 0.861eV
Energy per atom: Epot = 0.867eV  Ekin = 0.006eV (T= 47K)  Etot = 0.873eV
Energy per atom: Epot = 0.704eV  Ekin = 0.160eV (T=1241K)  Etot = 0.864eV
Energy per atom: Epot = 0.679eV  Ekin = 0.190eV (T=1468K)  Etot = 0.868eV
Energy per atom: Epot = 0.822eV  Ekin = 0.051eV (T=391K)  Etot = 0.872eV
Energy per atom: Epot = 0.782eV  Ekin = 0.088eV (T=677K)  Etot = 0.869eV
Energy per atom: Epot = 0.775eV  Ekin = 0.096eV (T=739K)  Etot = 0.871eV
Energy per atom: Epot = 0.751eV  Ekin = 0.119eV (T=920K)  Etot = 0.870eV
Energy per atom: Epot = 0.752eV  Ekin = 0.117eV (T=904K)  Etot = 0.869eV
Energy per atom: Epot = 0.764eV  Ekin = 0.107eV (T=825K)  Etot = 0.871eV
Energy per atom: Epot = 0.752eV  Ekin = 0.120eV (T=927K)  Etot = 0.871eV
Energy per atom: Epot = 0.747eV  Ekin = 0.123eV (T=953K)  Etot = 0.871eV
Energy per atom: Epot = 0.724eV  Ekin = 0.145eV (T=1119K)  Etot = 0.868eV
Energy per atom: Epot = 0.703eV  Ekin = 0.166eV (T=1281K)  Etot = 0.869eV
Energy per atom: Epot = 0.738eV  Ekin = 0.132eV (T=1022K)  Etot = 0.870eV
Energy per atom: Epot = 0.703eV  Ekin = 0.165eV (T=1280K)  Etot = 0.868eV
Energy per atom: Epot = 0.697eV  Ekin = 0.172eV (T=1331K)  Etot = 0.869eV
Energy per atom: Epot = 0.731eV  Ekin = 0.140eV (T=1083K)  Etot = 0.871eV
Energy per atom: Epot = 0.716eV  Ekin = 0.155eV (T=1199K)  Etot = 0.871eV
Energy per atom: Epot = 0.706eV  Ekin = 0.163eV (T=1257K)  Etot = 0.869eV
Energy per atom: Epot = 0.721eV  Ekin = 0.149eV (T=1151K)  Etot = 0.870eV
Energy per atom: Epot = 0.685eV  Ekin = 0.184eV (T=1424K)  Etot = 0.869eV
Energy per atom: Epot = 0.706eV  Ekin = 0.164eV (T=1271K)  Etot = 0.870eV
Energy per atom: Epot = 0.721eV  Ekin = 0.150eV (T=1161K)  Etot = 0.871eV
Energy per atom: Epot = 0.690eV  Ekin = 0.180eV (T=1396K)  Etot = 0.870eV
Energy per atom: Epot = 0.663eV  Ekin = 0.207eV (T=1601K)  Etot = 0.870eV
Energy per atom: Epot = 0.623eV  Ekin = 0.247eV (T=1907K)  Etot = 0.869eV
Energy per atom: Epot = 0.609eV  Ekin = 0.261eV (T=2022K)  Etot = 0.870eV
Energy per atom: Epot = 0.576eV  Ekin = 0.294eV (T=2271K)  Etot = 0.870eV
Energy per atom: Epot = 0.563eV  Ekin = 0.306eV (T=2370K)  Etot = 0.870eV
Energy per atom: Epot = 0.571eV  Ekin = 0.298eV (T=2307K)  Etot = 0.869eV
Energy per atom: Epot = 0.605eV  Ekin = 0.265eV (T=2050K)  Etot = 0.870eV
Energy per atom: Epot = 0.566eV  Ekin = 0.302eV (T=2338K)  Etot = 0.869eV
Energy per atom: Epot = 0.579eV  Ekin = 0.290eV (T=2242K)  Etot = 0.868eV
Energy per atom: Epot = 0.580eV  Ekin = 0.289eV (T=2234K)  Etot = 0.868eV
Energy per atom: Epot = 0.557eV  Ekin = 0.311eV (T=2406K)  Etot = 0.868eV
Energy per atom: Epot = 0.605eV  Ekin = 0.264eV (T=2046K)  Etot = 0.870eV
Energy per atom: Epot = 0.541eV  Ekin = 0.327eV (T=2530K)  Etot = 0.868eV
Energy per atom: Epot = 0.521eV  Ekin = 0.347eV (T=2686K)  Etot = 0.868eV
Energy per atom: Epot = 0.535eV  Ekin = 0.339eV (T=2622K)  Etot = 0.874eV
Energy per atom: Epot = 0.540eV  Ekin = 0.333eV (T=2574K)  Etot = 0.873eV
Energy per atom: Epot = 0.597eV  Ekin = 0.278eV (T=2151K)  Etot = 0.875eV
Energy per atom: Epot = 0.534eV  Ekin = 0.340eV (T=2631K)  Etot = 0.874eV
Energy per atom: Epot = 0.490eV  Ekin = 0.384eV (T=2968K)  Etot = 0.874eV
Energy per atom: Epot = 0.496eV  Ekin = 0.378eV (T=2928K)  Etot = 0.875eV
Energy per atom: Epot = 0.495eV  Ekin = 0.384eV (T=2970K)  Etot = 0.879eV
Energy per atom: Epot = 0.456eV  Ekin = 0.424eV (T=3277K)  Etot = 0.880eV
Energy per atom: Epot = 0.454eV  Ekin = 0.430eV (T=3325K)  Etot = 0.883eV
Energy per atom: Epot = 0.527eV  Ekin = 0.358eV (T=2767K)  Etot = 0.884eV
Energy per atom: Epot = 0.493eV  Ekin = 0.391eV (T=3027K)  Etot = 0.885eV
Energy per atom: Epot = 0.496eV  Ekin = 0.389eV (T=3011K)  Etot = 0.885eV
Energy per atom: Epot = 0.500eV  Ekin = 0.384eV (T=2971K)  Etot = 0.884eV
Energy per atom: Epot = 0.508eV  Ekin = 0.374eV (T=2895K)  Etot = 0.882eV
Energy per atom: Epot = 0.470eV  Ekin = 0.411eV (T=3176K)  Etot = 0.881eV
Energy per atom: Epot = 0.430eV  Ekin = 0.451eV (T=3489K)  Etot = 0.881eV
Energy per atom: Epot = 0.404eV  Ekin = 0.478eV (T=3696K)  Etot = 0.882eV
Energy per atom: Epot = 0.410eV  Ekin = 0.471eV (T=3645K)  Etot = 0.881eV
Energy per atom: Epot = 0.473eV  Ekin = 0.408eV (T=3160K)  Etot = 0.881eV
Energy per atom: Epot = 0.429eV  Ekin = 0.451eV (T=3490K)  Etot = 0.880eV
Energy per atom: Epot = 0.439eV  Ekin = 0.442eV (T=3423K)  Etot = 0.881eV
Energy per atom: Epot = 0.465eV  Ekin = 0.416eV (T=3216K)  Etot = 0.881eV
Energy per atom: Epot = 0.448eV  Ekin = 0.432eV (T=3346K)  Etot = 0.880eV
Energy per atom: Epot = 0.465eV  Ekin = 0.417eV (T=3222K)  Etot = 0.882eV
Energy per atom: Epot = 0.457eV  Ekin = 0.423eV (T=3274K)  Etot = 0.880eV
Energy per atom: Epot = 0.425eV  Ekin = 0.456eV (T=3528K)  Etot = 0.881eV
Energy per atom: Epot = 0.415eV  Ekin = 0.465eV (T=3598K)  Etot = 0.880eV
Energy per atom: Epot = 0.458eV  Ekin = 0.427eV (T=3301K)  Etot = 0.884eV
Energy per atom: Epot = 0.444eV  Ekin = 0.442eV (T=3423K)  Etot = 0.886eV
Energy per atom: Epot = 0.386eV  Ekin = 0.499eV (T=3858K)  Etot = 0.885eV
Energy per atom: Epot = 0.386eV  Ekin = 0.499eV (T=3864K)  Etot = 0.885eV
Energy per atom: Epot = 0.421eV  Ekin = 0.466eV (T=3609K)  Etot = 0.887eV
Energy per atom: Epot = 0.390eV  Ekin = 0.498eV (T=3854K)  Etot = 0.888eV
Energy per atom: Epot = 0.404eV  Ekin = 0.488eV (T=3777K)  Etot = 0.892eV
Energy per atom: Epot = 0.411eV  Ekin = 0.480eV (T=3714K)  Etot = 0.891eV
Energy per atom: Epot = 0.429eV  Ekin = 0.460eV (T=3558K)  Etot = 0.889eV
Energy per atom: Epot = 0.391eV  Ekin = 0.497eV (T=3848K)  Etot = 0.888eV
Energy per atom: Epot = 0.379eV  Ekin = 0.507eV (T=3924K)  Etot = 0.886eV
Energy per atom: Epot = 0.429eV  Ekin = 0.460eV (T=3555K)  Etot = 0.889eV
Energy per atom: Epot = 0.431eV  Ekin = 0.458eV (T=3545K)  Etot = 0.889eV
Energy per atom: Epot = 0.440eV  Ekin = 0.447eV (T=3458K)  Etot = 0.887eV
Energy per atom: Epot = 0.417eV  Ekin = 0.467eV (T=3610K)  Etot = 0.884eV
Energy per atom: Epot = 0.445eV  Ekin = 0.442eV (T=3422K)  Etot = 0.888eV
Energy per atom: Epot = 0.437eV  Ekin = 0.447eV (T=3461K)  Etot = 0.884eV
Energy per atom: Epot = 0.439eV  Ekin = 0.447eV (T=3460K)  Etot = 0.886eV
Energy per atom: Epot = 0.424eV  Ekin = 0.464eV (T=3587K)  Etot = 0.888eV
Energy per atom: Epot = 0.394eV  Ekin = 0.493eV (T=3817K)  Etot = 0.887eV
Energy per atom: Epot = 0.370eV  Ekin = 0.515eV (T=3984K)  Etot = 0.885eV
Energy per atom: Epot = 0.373eV  Ekin = 0.511eV (T=3951K)  Etot = 0.883eV
Energy per atom: Epot = 0.402eV  Ekin = 0.485eV (T=3748K)  Etot = 0.886eV
Energy per atom: Epot = 0.380eV  Ekin = 0.506eV (T=3912K)  Etot = 0.886eV
Energy per atom: Epot = 0.404eV  Ekin = 0.483eV (T=3736K)  Etot = 0.887eV
Energy per atom: Epot = 0.443eV  Ekin = 0.444eV (T=3436K)  Etot = 0.887eV
Energy per atom: Epot = 0.419eV  Ekin = 0.469eV (T=3626K)  Etot = 0.888eV
Energy per atom: Epot = 0.440eV  Ekin = 0.451eV (T=3489K)  Etot = 0.891eV
Energy per atom: Epot = 0.434eV  Ekin = 0.453eV (T=3505K)  Etot = 0.887eV
Energy per atom: Epot = 0.423eV  Ekin = 0.468eV (T=3619K)  Etot = 0.891eV
Energy per atom: Epot = 0.369eV  Ekin = 0.523eV (T=4047K)  Etot = 0.892eV
Energy per atom: Epot = 0.413eV  Ekin = 0.478eV (T=3701K)  Etot = 0.891eV
Energy per atom: Epot = 0.359eV  Ekin = 0.531eV (T=4108K)  Etot = 0.890eV
Energy per atom: Epot = 0.401eV  Ekin = 0.491eV (T=3799K)  Etot = 0.892eV
Energy per atom: Epot = 0.414eV  Ekin = 0.477eV (T=3688K)  Etot = 0.890eV
Energy per atom: Epot = 0.376eV  Ekin = 0.513eV (T=3969K)  Etot = 0.889eV
Energy per atom: Epot = 0.368eV  Ekin = 0.522eV (T=4035K)  Etot = 0.890eV
Energy per atom: Epot = 0.370eV  Ekin = 0.521eV (T=4029K)  Etot = 0.891eV
Energy per atom: Epot = 0.379eV  Ekin = 0.512eV (T=3959K)  Etot = 0.891eV
Energy per atom: Epot = 0.393eV  Ekin = 0.502eV (T=3881K)  Etot = 0.894eV
Energy per atom: Epot = 0.373eV  Ekin = 0.515eV (T=3984K)  Etot = 0.888eV
Energy per atom: Epot = 0.367eV  Ekin = 0.530eV (T=4100K)  Etot = 0.897eV
Energy per atom: Epot = 0.380eV  Ekin = 0.520eV (T=4026K)  Etot = 0.900eV
Energy per atom: Epot = 0.384eV  Ekin = 0.523eV (T=4048K)  Etot = 0.908eV
Energy per atom: Epot = 0.417eV  Ekin = 0.485eV (T=3751K)  Etot = 0.902eV
Energy per atom: Epot = 0.402eV  Ekin = 0.518eV (T=4005K)  Etot = 0.919eV
Energy per atom: Epot = 0.417eV  Ekin = 0.558eV (T=4313K)  Etot = 0.974eV
</pre></div></div>
</div>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[4]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="c1"># wrap and save traj in .xyz --- the .traj is a non human readable database file, xyz is much better</span>
<span class="n">out_traj</span> <span class="o">=</span> <span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">read</span><span class="p">(</span><span class="s1">&#39;dyn_emt.traj&#39;</span><span class="p">,</span> <span class="s1">&#39;:&#39;</span><span class="p">)</span>
<span class="k">for</span> <span class="n">at</span> <span class="ow">in</span> <span class="n">out_traj</span><span class="p">:</span>
    <span class="n">at</span><span class="o">.</span><span class="n">wrap</span><span class="p">()</span>
    <span class="k">if</span> <span class="s1">&#39;momenta&#39;</span> <span class="ow">in</span> <span class="n">at</span><span class="o">.</span><span class="n">arrays</span><span class="p">:</span> <span class="k">del</span> <span class="n">at</span><span class="o">.</span><span class="n">arrays</span><span class="p">[</span><span class="s1">&#39;momenta&#39;</span><span class="p">]</span>
<span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s1">&#39;dyn_emt.xyz&#39;</span><span class="p">,</span> <span class="n">out_traj</span><span class="p">,</span> <span class="s1">&#39;xyz&#39;</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
<div class="section" id="get-e0-for-H-and-O---energies-of-the-isolated-atoms">
<h2>get e0 for H and O - energies of the isolated atoms<a class="headerlink" href="#get-e0-for-H-and-O---energies-of-the-isolated-atoms" title="Permalink to this headline">¶</a></h2>
<p>This is the energy of the isolated atom, will be in the teach_sparse string in the following format: <code class="docutils literal notranslate"><span class="pre">e0={H:energy:O:energy}</span></code></p>
<div class="nbinput docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[5]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="n">isolated_H</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="s1">&#39;H&#39;</span><span class="p">,</span> <span class="n">calculator</span><span class="o">=</span><span class="n">EMT</span><span class="p">())</span>
<span class="n">isolated_O</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="s1">&#39;O&#39;</span><span class="p">,</span> <span class="n">calculator</span><span class="o">=</span><span class="n">EMT</span><span class="p">())</span>

<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;e0_H:&#39;</span><span class="p">,</span><span class="n">isolated_H</span><span class="o">.</span><span class="n">get_potential_energy</span><span class="p">())</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;e0_O:&#39;</span><span class="p">,</span><span class="n">isolated_O</span><span class="o">.</span><span class="n">get_potential_energy</span><span class="p">())</span>

<span class="c1"># this made the e0 string be the following: e0={H:3.21:O:4.6}</span>
</pre></div>
</div>
</div>
<div class="nboutput nblast docutils container">
<div class="prompt empty docutils container">
</div>
<div class="output_area docutils container">
<div class="highlight"><pre>
e0_H: 3.21
e0_O: 4.6
</pre></div></div>
</div>
</div>
</div>
<div class="section" id="separate-the-dataset-into-a-training-and-a-validation-set">
<h1>separate the dataset into a training and a validation set<a class="headerlink" href="#separate-the-dataset-into-a-training-and-a-validation-set" title="Permalink to this headline">¶</a></h1>
<p>As we have 120 frames from the 600fs MD, I will do it 60,60 with taking even and odd frames for the two</p>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[6]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s1">&#39;train.xyz&#39;</span><span class="p">,</span> <span class="n">out_traj</span><span class="p">[</span><span class="mi">0</span><span class="p">::</span><span class="mi">2</span><span class="p">])</span>
<span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s1">&#39;validate.xyz&#39;</span><span class="p">,</span> <span class="n">out_traj</span><span class="p">[</span><span class="mi">1</span><span class="p">::</span><span class="mi">2</span><span class="p">])</span>
</pre></div>
</div>
</div>
<div class="nbinput docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[7]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="n">out_traj</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">arrays</span><span class="o">.</span><span class="n">keys</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="nboutput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[7]:
</pre></div>
</div>
<div class="output_area highlight-none notranslate"><div class="highlight"><pre>
<span></span>dict_keys([&#39;numbers&#39;, &#39;positions&#39;])
</pre></div>
</div>
</div>
<div class="section" id="train-our-GAP-model-from-the-command-line">
<h2>train our GAP model from the command line<a class="headerlink" href="#train-our-GAP-model-from-the-command-line" title="Permalink to this headline">¶</a></h2>
<p>Will use a fit of 2b only, using the desciptor distance_2b.</p>
<p>Let’s understand how this works. The bash command takes named arguments separated by spaces.</p>
<ul class="simple">
<li><code class="docutils literal notranslate"><span class="pre">teach_sparse</span></code> the command which actually does the fit</li>
<li><code class="docutils literal notranslate"><span class="pre">e0={H:3.21:O:4.6}</span></code> the energies of the isolated atoms</li>
<li><code class="docutils literal notranslate"><span class="pre">energy_parameter_name=energy</span> <span class="pre">force_parameter_name=forces</span></code> names of the parameters</li>
<li><code class="docutils literal notranslate"><span class="pre">do_copy_at_file=F</span> <span class="pre">sparse_separate_file=T</span></code> just needed, don’t want to copy the training data and using separate files for the xml makes it faster</li>
<li><code class="docutils literal notranslate"><span class="pre">gp_file=GAP.xml</span></code> filename of the potential parameters, I have always used this name, because I had separate directories for the different trainings potentials</li>
<li><code class="docutils literal notranslate"><span class="pre">at_file=train.xyz</span></code> training file</li>
<li><code class="docutils literal notranslate"><span class="pre">default_sigma={0.008</span> <span class="pre">0.04</span> <span class="pre">0</span> <span class="pre">0}</span></code> sigma values to be used for energies, forces, stresses, hessians in order; this represents the accuracy of the data and the relative weight of them in the fit (more accurate –&gt; more significant in the fit)</li>
<li><code class="docutils literal notranslate"><span class="pre">gap={...}</span></code> the potential to be fit, separated by ‘:’</li>
</ul>
<p><strong>distance_2b</strong> - <code class="docutils literal notranslate"><span class="pre">cutoff=4.0</span></code> radial, practically the highest distance the descriptor takes into account - <code class="docutils literal notranslate"><span class="pre">covariance_type=ard_se</span></code> use gausses in the fit - <code class="docutils literal notranslate"><span class="pre">delta=0.5</span></code> what relative portion of the things shall be determined by this potential - <code class="docutils literal notranslate"><span class="pre">theta_uniform=1.0</span></code> width of the gaussians - <code class="docutils literal notranslate"><span class="pre">sparse_method=uniform</span></code> use uniform bins to choose the sparse points - <code class="docutils literal notranslate"><span class="pre">add_species=T</span></code> take the species into account, so it will generate more GAPs automatically (see the output) - <code class="docutils literal notranslate"><span class="pre">n_sparse=10</span></code>
number of sparse points</p>
</div>
<div class="section" id="notice,-that-the-script-is-running-in-parallel,-using-all-8-cores-of-the-current-machine">
<h2>notice, that the script is running in parallel, using all 8 cores of the current machine<a class="headerlink" href="#notice,-that-the-script-is-running-in-parallel,-using-all-8-cores-of-the-current-machine" title="Permalink to this headline">¶</a></h2>
<div class="nbinput docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[8]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="o">!</span> teach_sparse <span class="nv">e0</span><span class="o">={</span>H:3.21:O:4.6<span class="o">}</span> <span class="nv">energy_parameter_name</span><span class="o">=</span>energy <span class="nv">force_parameter_name</span><span class="o">=</span>forces <span class="nv">do_copy_at_file</span><span class="o">=</span>F <span class="nv">sparse_separate_file</span><span class="o">=</span>T <span class="nv">gp_file</span><span class="o">=</span>GAP.xml <span class="nv">at_file</span><span class="o">=</span>train.xyz <span class="nv">default_sigma</span><span class="o">={</span><span class="m">0</span>.008 <span class="m">0</span>.04 <span class="m">0</span> <span class="m">0</span><span class="o">}</span> <span class="nv">gap</span><span class="o">={</span>distance_2b <span class="nv">cutoff</span><span class="o">=</span><span class="m">4</span>.0 <span class="nv">covariance_type</span><span class="o">=</span>ard_se <span class="nv">delta</span><span class="o">=</span><span class="m">0</span>.5 <span class="nv">theta_uniform</span><span class="o">=</span><span class="m">1</span>.0 <span class="nv">sparse_method</span><span class="o">=</span>uniform <span class="nv">add_species</span><span class="o">=</span>T <span class="nv">n_sparse</span><span class="o">=</span><span class="m">10</span><span class="o">}</span>

</pre></div>
</div>
</div>
<div class="nboutput nblast docutils container">
<div class="prompt empty docutils container">
</div>
<div class="output_area docutils container">
<div class="highlight"><pre>
libAtoms::Hello World: 07/10/2018   18:10:52
libAtoms::Hello World: git version  https://github.com/libAtoms/QUIP.git,531330f-dirty
libAtoms::Hello World: QUIP_ARCH    linux_x86_64_gfortran_openmp
libAtoms::Hello World: compiled on  Jul  2 2018 at 21:44:13
libAtoms::Hello World: OpenMP parallelisation with 8 threads
WARNING: libAtoms::Hello World: environment variable OMP_STACKSIZE not set explicitly. The default value - system and compiler dependent - may be too small for some applications.
libAtoms::Hello World: Random Seed = 65452599
libAtoms::Hello World: global verbosity = 0

Calls to system_timer will do nothing by default


================================ Input parameters ==============================

at_file = train.xyz
gap = &#34;distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform add_species=T n_sparse=10&#34;
e0 = H:3.21:O:4.6
e0_offset = 0.0
do_e0_avg = T
default_sigma = &#34;0.008 0.04 0 0&#34;
sparse_jitter = 1.0e-10
hessian_delta = 1.0e-2
core_param_file = quip_params.xml
core_ip_args =
energy_parameter_name = energy
force_parameter_name = forces
virial_parameter_name = virial
hessian_parameter_name = hessian
config_type_parameter_name = config_type
sigma_parameter_name = sigma
config_type_sigma =
sigma_per_atom = T
do_copy_at_file = F
sparse_separate_file = T
sparse_use_actual_gpcov = F
gp_file = GAP.xml
verbosity = NORMAL
rnd_seed = -1
do_ip_timing = F
template_file = template.xyz

========================================  ======================================


============== Gaussian Approximation Potentials - Database fitting ============


Initial parsing of command line arguments finished.
Found 1 GAPs.
Descriptors have been parsed
XYZ file read
Old GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform add_species=T n_sparse=10}
New GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform               n_sparse=10 Z1=8 Z2=8}
New GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform               n_sparse=10 Z1=8 Z2=1}
New GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform               n_sparse=10 Z1=1 Z2=1}
Multispecies support added where requested
Number of target energies (property name: energy) found: 60
Number of target forces (property name: forces) found: 14580
Number of target virials (property name: virial) found: 0
Number of target Hessian eigenvalues (property name: hessian) found: 0
Cartesian coordinates transformed to descriptors
Started sparse covariance matrix calculation of coordinate 1

Finished sparse covariance matrix calculation of coordinate 1
TIMER: gpFull_covarianceMatrix_sparse_Coordinate1_sparse  done in .91605700000000012 cpu secs, .11785597354173660 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate1         done in .91605700000000012 cpu secs, .11797238700091839 wall clock secs.
Started sparse covariance matrix calculation of coordinate 2

Finished sparse covariance matrix calculation of coordinate 2
TIMER: gpFull_covarianceMatrix_sparse_Coordinate2_sparse  done in 3.7722349999999998 cpu secs, .47368377633392811 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate2         done in 3.7722349999999998 cpu secs, .47379869595170021 wall clock secs.
Started sparse covariance matrix calculation of coordinate 3

Finished sparse covariance matrix calculation of coordinate 3
TIMER: gpFull_covarianceMatrix_sparse_Coordinate3_sparse  done in 3.5482220000000000 cpu secs, .44678380154073238 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate3         done in 3.5522219999999995 cpu secs, .44689681194722652 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_LinearAlgebra       done in .56005000000000749E-001 cpu secs, .13479800894856453E-001 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_FunctionValues      done in .00000000000000000E+000 cpu secs, .14016032218933105E-003 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse                     done in 8.3005189999999995 cpu secs, 1.0586479119956493 wall clock secs.
TIMER: GP sparsify                                        done in 9.8286140000000000 cpu secs, 1.4589964970946312 wall clock secs.

libAtoms::Finalise: 07/10/2018   18:10:55
libAtoms::Finalise: Bye-Bye!
</pre></div></div>
</div>
</div>
<div class="section" id="use-the-potential-with-QUIP-on-trani.xyz-and-validate.xyz">
<h2>use the potential with QUIP on trani.xyz and validate.xyz<a class="headerlink" href="#use-the-potential-with-QUIP-on-trani.xyz-and-validate.xyz" title="Permalink to this headline">¶</a></h2>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[9]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="c1"># calculate train.xyz</span>

<span class="o">!</span> quip <span class="nv">E</span><span class="o">=</span>T <span class="nv">F</span><span class="o">=</span>T <span class="nv">atoms_filename</span><span class="o">=</span>train.xyz <span class="nv">param_filename</span><span class="o">=</span>GAP.xml <span class="p">|</span> grep AT <span class="p">|</span> sed <span class="s1">&#39;s/AT//&#39;</span> &gt;&gt; quip_train.xyz
<span class="o">!</span> quip <span class="nv">E</span><span class="o">=</span>T <span class="nv">F</span><span class="o">=</span>T <span class="nv">atoms_filename</span><span class="o">=</span>validate.xyz <span class="nv">param_filename</span><span class="o">=</span>GAP.xml <span class="p">|</span> grep AT <span class="p">|</span> sed <span class="s1">&#39;s/AT//&#39;</span> &gt;&gt; quip_validate.xyz
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="make-simple-plots-of-the-energies-and-forces-on-the-EMT-and-GAP-datas">
<h1>make simple plots of the energies and forces on the EMT and GAP datas<a class="headerlink" href="#make-simple-plots-of-the-energies-and-forces-on-the-EMT-and-GAP-datas" title="Permalink to this headline">¶</a></h1>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[16]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="k">def</span> <span class="nf">energy_plot</span><span class="p">(</span><span class="n">in_file</span><span class="p">,</span> <span class="n">out_file</span><span class="p">,</span> <span class="n">ax</span><span class="p">,</span> <span class="n">title</span><span class="o">=</span><span class="s1">&#39;Plot of energy&#39;</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot; Plots the distribution of energy per atom on the output vs the input&quot;&quot;&quot;</span>
    <span class="c1"># read files</span>
    <span class="n">in_atoms</span> <span class="o">=</span> <span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">read</span><span class="p">(</span><span class="n">in_file</span><span class="p">,</span> <span class="s1">&#39;:&#39;</span><span class="p">)</span>
    <span class="n">out_atoms</span> <span class="o">=</span> <span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">read</span><span class="p">(</span><span class="n">out_file</span><span class="p">,</span> <span class="s1">&#39;:&#39;</span><span class="p">)</span>
    <span class="c1"># list energies</span>
    <span class="n">ener_in</span> <span class="o">=</span> <span class="p">[</span><span class="n">at</span><span class="o">.</span><span class="n">get_potential_energy</span><span class="p">()</span> <span class="o">/</span> <span class="nb">len</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">get_chemical_symbols</span><span class="p">())</span> <span class="k">for</span> <span class="n">at</span> <span class="ow">in</span> <span class="n">in_atoms</span><span class="p">]</span>
    <span class="n">ener_out</span> <span class="o">=</span> <span class="p">[</span><span class="n">at</span><span class="o">.</span><span class="n">get_potential_energy</span><span class="p">()</span> <span class="o">/</span> <span class="nb">len</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">get_chemical_symbols</span><span class="p">())</span> <span class="k">for</span> <span class="n">at</span> <span class="ow">in</span> <span class="n">out_atoms</span><span class="p">]</span>
    <span class="c1"># scatter plot of the data</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">ener_in</span><span class="p">,</span> <span class="n">ener_out</span><span class="p">)</span>
    <span class="c1"># get the appropriate limits for the plot</span>
    <span class="n">for_limits</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">ener_in</span> <span class="o">+</span><span class="n">ener_out</span><span class="p">)</span>
    <span class="n">elim</span> <span class="o">=</span> <span class="p">(</span><span class="n">for_limits</span><span class="o">.</span><span class="n">min</span><span class="p">()</span> <span class="o">-</span> <span class="mf">0.05</span><span class="p">,</span> <span class="n">for_limits</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">+</span> <span class="mf">0.05</span><span class="p">)</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">(</span><span class="n">elim</span><span class="p">)</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">(</span><span class="n">elim</span><span class="p">)</span>
    <span class="c1"># add line of slope 1 for refrence</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">elim</span><span class="p">,</span> <span class="n">elim</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="s1">&#39;k&#39;</span><span class="p">)</span>
    <span class="c1"># set labels</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;energy by GAP / eV&#39;</span><span class="p">)</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;energy by EMT / eV&#39;</span><span class="p">)</span>
    <span class="c1">#set title</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_title</span><span class="p">(</span><span class="n">title</span><span class="p">)</span>
    <span class="c1"># add text about RMSE</span>
    <span class="n">_rms</span> <span class="o">=</span> <span class="n">rms_dict</span><span class="p">(</span><span class="n">ener_in</span><span class="p">,</span> <span class="n">ener_out</span><span class="p">)</span>
    <span class="n">rmse_text</span> <span class="o">=</span> <span class="s1">&#39;RMSE:</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">round</span><span class="p">(</span><span class="n">_rms</span><span class="p">[</span><span class="s1">&#39;rmse&#39;</span><span class="p">],</span> <span class="mi">3</span><span class="p">))</span> <span class="o">+</span> <span class="s1">&#39; +- &#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">round</span><span class="p">(</span><span class="n">_rms</span><span class="p">[</span><span class="s1">&#39;std&#39;</span><span class="p">],</span> <span class="mi">3</span><span class="p">))</span> <span class="o">+</span> <span class="s1">&#39;eV/atom&#39;</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="mf">0.9</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">,</span> <span class="n">rmse_text</span><span class="p">,</span> <span class="n">transform</span><span class="o">=</span><span class="n">ax</span><span class="o">.</span><span class="n">transAxes</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="s1">&#39;large&#39;</span><span class="p">,</span> <span class="n">horizontalalignment</span><span class="o">=</span><span class="s1">&#39;right&#39;</span><span class="p">,</span>
            <span class="n">verticalalignment</span><span class="o">=</span><span class="s1">&#39;bottom&#39;</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">force_plot</span><span class="p">(</span><span class="n">in_file</span><span class="p">,</span> <span class="n">out_file</span><span class="p">,</span> <span class="n">ax</span><span class="p">,</span> <span class="n">symbol</span><span class="o">=</span><span class="s1">&#39;HO&#39;</span><span class="p">,</span> <span class="n">title</span><span class="o">=</span><span class="s1">&#39;Plot of force&#39;</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot; Plots the distribution of firce components per atom on the output vs the input</span>
<span class="sd">        only plots for the given atom type(s)&quot;&quot;&quot;</span>

    <span class="n">in_atoms</span> <span class="o">=</span> <span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">read</span><span class="p">(</span><span class="n">in_file</span><span class="p">,</span> <span class="s1">&#39;:&#39;</span><span class="p">)</span>
    <span class="n">out_atoms</span> <span class="o">=</span> <span class="n">ase</span><span class="o">.</span><span class="n">io</span><span class="o">.</span><span class="n">read</span><span class="p">(</span><span class="n">out_file</span><span class="p">,</span> <span class="s1">&#39;:&#39;</span><span class="p">)</span>

    <span class="c1"># extract data for only one species</span>
    <span class="n">in_force</span><span class="p">,</span> <span class="n">out_force</span> <span class="o">=</span> <span class="p">[],</span> <span class="p">[]</span>
    <span class="k">for</span> <span class="n">at_in</span><span class="p">,</span> <span class="n">at_out</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="n">in_atoms</span><span class="p">,</span> <span class="n">out_atoms</span><span class="p">):</span>
        <span class="c1"># get the symbols</span>
        <span class="n">sym_all</span> <span class="o">=</span> <span class="n">at_in</span><span class="o">.</span><span class="n">get_chemical_symbols</span><span class="p">()</span>
        <span class="c1"># add force for each atom</span>
        <span class="k">for</span> <span class="n">j</span><span class="p">,</span> <span class="n">sym</span> <span class="ow">in</span> <span class="nb">enumerate</span><span class="p">(</span><span class="n">sym_all</span><span class="p">):</span>
            <span class="k">if</span> <span class="n">sym</span> <span class="ow">in</span> <span class="n">symbol</span><span class="p">:</span>
                <span class="n">in_force</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">at_in</span><span class="o">.</span><span class="n">get_forces</span><span class="p">()[</span><span class="n">j</span><span class="p">])</span>
                <span class="c1">#out_force.append(at_out.get_forces()[j]) \</span>
                <span class="n">out_force</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">at_out</span><span class="o">.</span><span class="n">arrays</span><span class="p">[</span><span class="s1">&#39;force&#39;</span><span class="p">][</span><span class="n">j</span><span class="p">])</span> <span class="c1"># because QUIP and ASE use different names</span>
    <span class="c1"># convert to np arrays, much easier to work with</span>
    <span class="c1">#in_force = np.array(in_force)</span>
    <span class="c1">#out_force = np.array(out_force)</span>
    <span class="c1"># scatter plot of the data</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">scatter</span><span class="p">(</span><span class="n">in_force</span><span class="p">,</span> <span class="n">out_force</span><span class="p">)</span>
    <span class="c1"># get the appropriate limits for the plot</span>
    <span class="n">for_limits</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">in_force</span> <span class="o">+</span> <span class="n">out_force</span><span class="p">)</span>
    <span class="n">flim</span> <span class="o">=</span> <span class="p">(</span><span class="n">for_limits</span><span class="o">.</span><span class="n">min</span><span class="p">()</span> <span class="o">-</span> <span class="mi">1</span><span class="p">,</span> <span class="n">for_limits</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">+</span> <span class="mi">1</span><span class="p">)</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_xlim</span><span class="p">(</span><span class="n">flim</span><span class="p">)</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_ylim</span><span class="p">(</span><span class="n">flim</span><span class="p">)</span>
    <span class="c1"># add line of</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">flim</span><span class="p">,</span> <span class="n">flim</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="s1">&#39;k&#39;</span><span class="p">)</span>
    <span class="c1"># set labels</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_ylabel</span><span class="p">(</span><span class="s1">&#39;force by GAP / (eV/Å)&#39;</span><span class="p">)</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_xlabel</span><span class="p">(</span><span class="s1">&#39;force by EMT / (eV/Å)&#39;</span><span class="p">)</span>
    <span class="c1">#set title</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">set_title</span><span class="p">(</span><span class="n">title</span><span class="p">)</span>
    <span class="c1"># add text about RMSE</span>
    <span class="n">_rms</span> <span class="o">=</span> <span class="n">rms_dict</span><span class="p">(</span><span class="n">in_force</span><span class="p">,</span> <span class="n">out_force</span><span class="p">)</span>
    <span class="n">rmse_text</span> <span class="o">=</span> <span class="s1">&#39;RMSE:</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">round</span><span class="p">(</span><span class="n">_rms</span><span class="p">[</span><span class="s1">&#39;rmse&#39;</span><span class="p">],</span> <span class="mi">3</span><span class="p">))</span> <span class="o">+</span> <span class="s1">&#39; +- &#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">round</span><span class="p">(</span><span class="n">_rms</span><span class="p">[</span><span class="s1">&#39;std&#39;</span><span class="p">],</span> <span class="mi">3</span><span class="p">))</span> <span class="o">+</span> <span class="s1">&#39;eV/Å&#39;</span>
    <span class="n">ax</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="mf">0.9</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">,</span> <span class="n">rmse_text</span><span class="p">,</span> <span class="n">transform</span><span class="o">=</span><span class="n">ax</span><span class="o">.</span><span class="n">transAxes</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="s1">&#39;large&#39;</span><span class="p">,</span> <span class="n">horizontalalignment</span><span class="o">=</span><span class="s1">&#39;right&#39;</span><span class="p">,</span>
            <span class="n">verticalalignment</span><span class="o">=</span><span class="s1">&#39;bottom&#39;</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="nbinput docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[25]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="n">fig</span><span class="p">,</span> <span class="n">ax_list</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="n">nrows</span><span class="o">=</span><span class="mi">3</span><span class="p">,</span> <span class="n">ncols</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">gridspec_kw</span><span class="o">=</span><span class="p">{</span><span class="s1">&#39;hspace&#39;</span><span class="p">:</span> <span class="mf">0.3</span><span class="p">})</span>
<span class="n">fig</span><span class="o">.</span><span class="n">set_size_inches</span><span class="p">(</span><span class="mi">15</span><span class="p">,</span> <span class="mi">20</span><span class="p">)</span>
<span class="n">ax_list</span> <span class="o">=</span> <span class="n">ax_list</span><span class="o">.</span><span class="n">flat</span><span class="p">[:]</span>

<span class="n">energy_plot</span><span class="p">(</span><span class="s1">&#39;train.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_train.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="s1">&#39;Energy on training data&#39;</span><span class="p">)</span>
<span class="n">energy_plot</span><span class="p">(</span><span class="s1">&#39;validate.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_validate.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="s1">&#39;Energy on validation data&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;train.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_train.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="s1">&#39;H&#39;</span><span class="p">,</span> <span class="s1">&#39;Force on training data - H&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;train.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_train.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span> <span class="s1">&#39;O&#39;</span><span class="p">,</span> <span class="s1">&#39;Force on training data - O&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;validate.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_validate.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span> <span class="s1">&#39;H&#39;</span><span class="p">,</span> <span class="s1">&#39;Force on validation data - H&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;validate.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_validate.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span> <span class="s1">&#39;O&#39;</span><span class="p">,</span>  <span class="s1">&#39;Force on validation data - O&#39;</span><span class="p">)</span>

<span class="c1"># if you wanted to have the same limits on the firce plots</span>
<span class="c1">#for ax in ax_list[2:]:</span>
<span class="c1">#    flim = (-20, 20)</span>
<span class="c1">#    ax.set_xlim(flim)</span>
<span class="c1">#    ax.set_ylim(flim)</span>
</pre></div>
</div>
</div>
<div class="nboutput nblast docutils container">
<div class="prompt empty docutils container">
</div>
<div class="output_area docutils container">
<img alt="../_images/Tutorials_gap_fitting_tutorial_17_0.png" src="../_images/Tutorials_gap_fitting_tutorial_17_0.png" />
</div>
</div>
<div class="section" id="train-our-GAP_3b-model-from-the-command-line">
<h2>train our GAP_3b model from the command line<a class="headerlink" href="#train-our-GAP_3b-model-from-the-command-line" title="Permalink to this headline">¶</a></h2>
<p>Let’s add three ody terms to the fit, which will hopefully improve it. We will be using the desciprtors distance_2b and angle_3b.</p>
<p><strong>angle_3b</strong> - <code class="docutils literal notranslate"><span class="pre">theta_fac=0.5</span></code> this takes the input data and determines the width from that; useful here, because the dimensions of the descriptor are different - <code class="docutils literal notranslate"><span class="pre">n_sparse=50</span></code> higher dimensional space, more sparse points</p>
</div>
<div class="section" id="both-training-and-quip-takes-significantly-more-time-than-the-last-one!!!">
<h2>both training and quip takes significantly more time than the last one!!!<a class="headerlink" href="#both-training-and-quip-takes-significantly-more-time-than-the-last-one!!!" title="Permalink to this headline">¶</a></h2>
<div class="nbinput docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[12]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="o">!</span> teach_sparse <span class="nv">e0</span><span class="o">={</span>H:3.21:O:4.6<span class="o">}</span> <span class="nv">energy_parameter_name</span><span class="o">=</span>energy <span class="nv">force_parameter_name</span><span class="o">=</span>forces <span class="nv">do_copy_at_file</span><span class="o">=</span>F <span class="nv">sparse_separate_file</span><span class="o">=</span>T <span class="nv">gp_file</span><span class="o">=</span>GAP_3b.xml <span class="nv">at_file</span><span class="o">=</span>train.xyz <span class="nv">default_sigma</span><span class="o">={</span><span class="m">0</span>.008 <span class="m">0</span>.04 <span class="m">0</span> <span class="m">0</span><span class="o">}</span> <span class="nv">gap</span><span class="o">={</span>distance_2b <span class="nv">cutoff</span><span class="o">=</span><span class="m">4</span>.0 <span class="nv">covariance_type</span><span class="o">=</span>ard_se <span class="nv">delta</span><span class="o">=</span><span class="m">0</span>.5 <span class="nv">theta_uniform</span><span class="o">=</span><span class="m">1</span>.0 <span class="nv">sparse_method</span><span class="o">=</span>uniform <span class="nv">add_species</span><span class="o">=</span>T <span class="nv">n_sparse</span><span class="o">=</span><span class="m">10</span> : angle_3b <span class="nv">cutoff</span><span class="o">=</span><span class="m">3</span>.5 <span class="nv">covariance_type</span><span class="o">=</span>ard_se <span class="nv">delta</span><span class="o">=</span><span class="m">0</span>.5 <span class="nv">theta_fac</span><span class="o">=</span><span class="m">0</span>.5 <span class="nv">add_species</span><span class="o">=</span>T <span class="nv">n_sparse</span><span class="o">=</span><span class="m">30</span> <span class="nv">sparse_method</span><span class="o">=</span>uniform<span class="o">}</span>

</pre></div>
</div>
</div>
<div class="nboutput nblast docutils container">
<div class="prompt empty docutils container">
</div>
<div class="output_area docutils container">
<div class="highlight"><pre>
libAtoms::Hello World: 07/10/2018   18:11:09
libAtoms::Hello World: git version  https://github.com/libAtoms/QUIP.git,531330f-dirty
libAtoms::Hello World: QUIP_ARCH    linux_x86_64_gfortran_openmp
libAtoms::Hello World: compiled on  Jul  2 2018 at 21:44:13
libAtoms::Hello World: OpenMP parallelisation with 8 threads
WARNING: libAtoms::Hello World: environment variable OMP_STACKSIZE not set explicitly. The default value - system and compiler dependent - may be too small for some applications.
libAtoms::Hello World: Random Seed = 65469770
libAtoms::Hello World: global verbosity = 0

Calls to system_timer will do nothing by default


================================ Input parameters ==============================

at_file = train.xyz
gap = &#34;distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform add_species=T n_sparse=10 : angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5 add_species=T n_sparse=30 sparse_method=uniform&#34;
e0 = H:3.21:O:4.6
e0_offset = 0.0
do_e0_avg = T
default_sigma = &#34;0.008 0.04 0 0&#34;
sparse_jitter = 1.0e-10
hessian_delta = 1.0e-2
core_param_file = quip_params.xml
core_ip_args =
energy_parameter_name = energy
force_parameter_name = forces
virial_parameter_name = virial
hessian_parameter_name = hessian
config_type_parameter_name = config_type
sigma_parameter_name = sigma
config_type_sigma =
sigma_per_atom = T
do_copy_at_file = F
sparse_separate_file = T
sparse_use_actual_gpcov = F
gp_file = GAP_3b.xml
verbosity = NORMAL
rnd_seed = -1
do_ip_timing = F
template_file = template.xyz

========================================  ======================================


============== Gaussian Approximation Potentials - Database fitting ============


Initial parsing of command line arguments finished.
Found 2 GAPs.
Descriptors have been parsed
XYZ file read
Old GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform add_species=T n_sparse=10}
New GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform               n_sparse=10 Z1=8 Z2=8}
New GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform               n_sparse=10 Z1=8 Z2=1}
New GAP: {distance_2b cutoff=4.0 covariance_type=ard_se delta=0.5 theta_uniform=1.0 sparse_method=uniform               n_sparse=10 Z1=1 Z2=1}
Old GAP: { angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5 add_species=T n_sparse=30 sparse_method=uniform}
New GAP: { angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5               n_sparse=30 sparse_method=uniform Z=8 Z1=8 Z2=8}
New GAP: { angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5               n_sparse=30 sparse_method=uniform Z=8 Z1=8 Z2=1}
New GAP: { angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5               n_sparse=30 sparse_method=uniform Z=8 Z1=1 Z2=1}
New GAP: { angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5               n_sparse=30 sparse_method=uniform Z=1 Z1=8 Z2=8}
New GAP: { angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5               n_sparse=30 sparse_method=uniform Z=1 Z1=8 Z2=1}
New GAP: { angle_3b cutoff=3.5 covariance_type=ard_se delta=0.5 theta_fac=0.5               n_sparse=30 sparse_method=uniform Z=1 Z1=1 Z2=1}
Multispecies support added where requested
Number of target energies (property name: energy) found: 60
Number of target forces (property name: forces) found: 14580
Number of target virials (property name: virial) found: 0
Number of target Hessian eigenvalues (property name: hessian) found: 0
Cartesian coordinates transformed to descriptors
Started sparse covariance matrix calculation of coordinate 1

Finished sparse covariance matrix calculation of coordinate 1
TIMER: gpFull_covarianceMatrix_sparse_Coordinate1_sparse  done in .89205500000000093 cpu secs, .11343244090676308 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate1         done in .89205500000000093 cpu secs, .11355240270495415 wall clock secs.
Started sparse covariance matrix calculation of coordinate 2

Finished sparse covariance matrix calculation of coordinate 2
TIMER: gpFull_covarianceMatrix_sparse_Coordinate2_sparse  done in 3.6442280000000018 cpu secs, .45859701186418533 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate2         done in 3.6442280000000018 cpu secs, .45871092379093170 wall clock secs.
Started sparse covariance matrix calculation of coordinate 3

Finished sparse covariance matrix calculation of coordinate 3
TIMER: gpFull_covarianceMatrix_sparse_Coordinate3_sparse  done in 3.4482159999999986 cpu secs, .43381043337285519 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate3         done in 3.4482159999999986 cpu secs, .43392318300902843 wall clock secs.
Started sparse covariance matrix calculation of coordinate 4

Finished sparse covariance matrix calculation of coordinate 4
TIMER: gpFull_covarianceMatrix_sparse_Coordinate4_sparse  done in 9.9606220000000008 cpu secs, 1.2572825513780117 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate4         done in 9.9606220000000008 cpu secs, 1.2574077267199755 wall clock secs.
Started sparse covariance matrix calculation of coordinate 5

Finished sparse covariance matrix calculation of coordinate 5
TIMER: gpFull_covarianceMatrix_sparse_Coordinate5_sparse  done in 50.367147000000003 cpu secs, 6.3618728630244732 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate5         done in 50.367147000000003 cpu secs, 6.3619935847818851 wall clock secs.
Started sparse covariance matrix calculation of coordinate 6

Finished sparse covariance matrix calculation of coordinate 6
TIMER: gpFull_covarianceMatrix_sparse_Coordinate6_sparse  done in 55.067442000000000 cpu secs, 6.9305843152105808 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate6         done in 55.067442000000000 cpu secs, 6.9306991323828697 wall clock secs.
Started sparse covariance matrix calculation of coordinate 7

Finished sparse covariance matrix calculation of coordinate 7
TIMER: gpFull_covarianceMatrix_sparse_Coordinate7_sparse  done in 25.993623999999983 cpu secs, 3.2758836857974529 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate7         done in 25.993623999999983 cpu secs, 3.2760034948587418 wall clock secs.
Started sparse covariance matrix calculation of coordinate 8

Finished sparse covariance matrix calculation of coordinate 8
TIMER: gpFull_covarianceMatrix_sparse_Coordinate8_sparse  done in 112.66704200000001 cpu secs, 14.186459138989449 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate8         done in 112.67104200000003 cpu secs, 14.186579804867506 wall clock secs.
Started sparse covariance matrix calculation of coordinate 9

Finished sparse covariance matrix calculation of coordinate 9
TIMER: gpFull_covarianceMatrix_sparse_Coordinate9_sparse  done in 103.17444799999998 cpu secs, 12.990166073665023 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_Coordinate9         done in 103.17444799999998 cpu secs, 12.990285659208894 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_LinearAlgebra       done in .44402700000000550 cpu secs, .12523609958589077 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse_FunctionValues      done in .00000000000000000E+000 cpu secs, .13986043632030487E-003 wall clock secs.
TIMER: gpFull_covarianceMatrix_sparse                     done in 365.71485500000000 cpu secs, 46.178485954180360 wall clock secs.
TIMER: GP sparsify                                        done in 368.54703300000000 cpu secs, 47.065343659371138 wall clock secs.

libAtoms::Finalise: 07/10/2018   18:12:12
libAtoms::Finalise: Bye-Bye!
</pre></div></div>
</div>
</div>
<div class="section" id="use-the-potential-with-QUIP-on-trani.xyz-and-validate.xyz">
<h2>use the potential with QUIP on trani.xyz and validate.xyz<a class="headerlink" href="#use-the-potential-with-QUIP-on-trani.xyz-and-validate.xyz" title="Permalink to this headline">¶</a></h2>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[13]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="c1"># calculate train.xyz</span>

<span class="o">!</span> quip <span class="nv">E</span><span class="o">=</span>T <span class="nv">F</span><span class="o">=</span>T <span class="nv">atoms_filename</span><span class="o">=</span>train.xyz <span class="nv">param_filename</span><span class="o">=</span>GAP_3b.xml <span class="p">|</span> grep AT <span class="p">|</span> sed <span class="s1">&#39;s/AT//&#39;</span> &gt;&gt; quip_3b_train.xyz
<span class="o">!</span> quip <span class="nv">E</span><span class="o">=</span>T <span class="nv">F</span><span class="o">=</span>T <span class="nv">atoms_filename</span><span class="o">=</span>validate.xyz <span class="nv">param_filename</span><span class="o">=</span>GAP_3b.xml <span class="p">|</span> grep AT <span class="p">|</span> sed <span class="s1">&#39;s/AT//&#39;</span> &gt;&gt; quip_3b_validate.xyz
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="look-at-the-outputs---clear-improvement">
<h1>look at the outputs - clear improvement<a class="headerlink" href="#look-at-the-outputs---clear-improvement" title="Permalink to this headline">¶</a></h1>
<div class="nbinput docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[26]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span><span class="n">fig</span><span class="p">,</span> <span class="n">ax_list</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">subplots</span><span class="p">(</span><span class="n">nrows</span><span class="o">=</span><span class="mi">3</span><span class="p">,</span> <span class="n">ncols</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">gridspec_kw</span><span class="o">=</span><span class="p">{</span><span class="s1">&#39;hspace&#39;</span><span class="p">:</span> <span class="mf">0.3</span><span class="p">})</span>
<span class="n">fig</span><span class="o">.</span><span class="n">set_size_inches</span><span class="p">(</span><span class="mi">15</span><span class="p">,</span> <span class="mi">20</span><span class="p">)</span>
<span class="n">ax_list</span> <span class="o">=</span> <span class="n">ax_list</span><span class="o">.</span><span class="n">flat</span><span class="p">[:]</span>

<span class="n">energy_plot</span><span class="p">(</span><span class="s1">&#39;train.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_3b_train.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="s1">&#39;Energy on training data&#39;</span><span class="p">)</span>
<span class="n">energy_plot</span><span class="p">(</span><span class="s1">&#39;validate.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_3b_validate.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="s1">&#39;Energy on validation data&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;train.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_3b_train.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="s1">&#39;H&#39;</span><span class="p">,</span> <span class="s1">&#39;Force on training data - H&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;train.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_3b_train.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span> <span class="s1">&#39;O&#39;</span><span class="p">,</span> <span class="s1">&#39;Force on training data - O&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;validate.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_3b_validate.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span> <span class="s1">&#39;H&#39;</span><span class="p">,</span> <span class="s1">&#39;Force on validation data - H&#39;</span><span class="p">)</span>
<span class="n">force_plot</span><span class="p">(</span><span class="s1">&#39;validate.xyz&#39;</span><span class="p">,</span> <span class="s1">&#39;quip_3b_validate.xyz&#39;</span><span class="p">,</span> <span class="n">ax_list</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span> <span class="s1">&#39;O&#39;</span><span class="p">,</span>  <span class="s1">&#39;Force on validation data - O&#39;</span><span class="p">)</span>

<span class="c1"># if you wanted to have the same limits on the firce plots</span>
<span class="c1">#for ax in ax_list[2:]:</span>
<span class="c1">#    flim = (-20, 20)</span>
<span class="c1">#    ax.set_xlim(flim)</span>
<span class="c1">#    ax.set_ylim(flim)</span>
</pre></div>
</div>
</div>
<div class="nboutput nblast docutils container">
<div class="prompt empty docutils container">
</div>
<div class="output_area docutils container">
<img alt="../_images/Tutorials_gap_fitting_tutorial_23_0.png" src="../_images/Tutorials_gap_fitting_tutorial_23_0.png" />
</div>
</div>
<div class="nbinput nblast docutils container">
<div class="prompt highlight-none notranslate"><div class="highlight"><pre><span></span>[ ]:
</pre></div>
</div>
<div class="input_area highlight-ipython3 notranslate"><div class="highlight"><pre>
<span></span>
</pre></div>
</div>
</div>
</div>


           </div>
           
          </div>
          <footer>
  

  <hr/>

  <div role="contentinfo">
    <p>
        &copy; Copyright 2008-2016, James Kermode

    </p>
  </div>
  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>. 

</footer>

        </div>
      </div>

    </section>

  </div>
  


  <script type="text/javascript">
      jQuery(function () {
          SphinxRtdTheme.Navigation.enable(true);
      });
  </script>

  
  
    
   

</body>
</html>